package com.numericalmethod.algoquant.model.ralph2009;

import org.apache.commons.math3.analysis.function.Power;
import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.ArrayRealVector;
import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.linear.RealVector;
import org.apache.commons.math3.stat.regression.OLSMultipleLinearRegression;

/**
 * Helper class for the power series regression described in the equation (19) and (20) in the paper
 * 
 * A Simulation Approach to Dynamic Portfolio Choice with an Application to Learning About Return Predictability
 * (Brandt 2005)
 *  
 * @author Paul/Clement/Stephen
 *
 */
public class Brandt2005Regression {
	
	/**
	 * Implementation of the estimator using the power series regression described in Brandt 2005 equation (19) and (20), 
	 * in the form of double array values instead of Vectors
	 */
	public static double[] estimateYValues(double[] z1, double[] z2, double[] y){
		return estimateYValues(new ArrayRealVector(z1), new ArrayRealVector(z2), new ArrayRealVector(y)).toArray();
	}
	
	/**
	 * Implementation of the estimator using the power series regression 
	 * the following formula described in Brandt 2005 equation (19) and (20). <br/>
	 * 
	 * <br/><b>
	 * 
	 * y = ( 1   z1  z2   z1^2  z2^2 ) * transpose (theta0 theta1 theta2 theta3 theta4)  
	 * 
	 * </b><br/>
	 * 
	 * @param vecZ1 The first independent variable observations for regression (z1)
	 * @param vecZ2 The second independent variable observations for regression (z2)
	 * @param vecY The dependent variable observations for regression (y)
	 * @return The estimated y using the regression result 
	 */
	public static RealVector estimateYValues(RealVector vecZ1, RealVector vecZ2, RealVector vecY){
		OLSMultipleLinearRegression regression = new OLSMultipleLinearRegression();
		RealVector vecZ1Sq = vecZ1.map(new Power(2));
		RealVector vecZ2Sq = vecZ2.map(new Power(2));
		RealVector unitVec = new ArrayRealVector(vecZ1.getDimension());
		unitVec.set(1);
	
		RealMatrix zMatrix = new Array2DRowRealMatrix(vecZ1.getDimension(),5);
		zMatrix.setColumnVector(0, unitVec); 
		//unary column is inserted internally
		zMatrix.setColumnVector(1, vecZ1);
		zMatrix.setColumnVector(2, vecZ2);
		zMatrix.setColumnVector(3, vecZ1Sq);
		zMatrix.setColumnVector(4, vecZ2Sq);
	
		RealMatrix subMatrixZ = zMatrix.getSubMatrix(0,zMatrix.getRowDimension()-1,1,zMatrix.getColumnDimension()-1);
		regression.newSampleData(vecY.toArray(), subMatrixZ.getData());
		double[] theta = regression.estimateRegressionParameters();
		
		RealMatrix thetaMatrix =  new Array2DRowRealMatrix(theta);

		RealVector estimatedYs = zMatrix.multiply(thetaMatrix).getColumnVector(0);
		
		return estimatedYs;
	}
	
	
	

}
